//File: geometry.hpp
//Date: Sun Oct 15 2021
//Author: LiPeng(m201572398@alumni.hust.edu.cn)

#pragma once

#include <vector>
#include <limits>
#include <algorithm>
#include <cmath>
#include "utils.hpp"
#include "debugutils.hpp"

class Geometry {
    public:
        int w, h;

        Geometry(int m_w=0, int m_h=0): w(m_w), h(m_h) 
        {}

        int area() const
        {
            return w * h;
        }

        real_t ratio() const
        {
            return (real_t) std::max(w, h) / std::min(w, h);
        }

        bool contain(int x, int y)
        {return(x >= 0) && (x < w) && (y >= 0) && (y < h);}        
};

template<typename T>
class Vector {
    public:
        T x = 0, y = 0, z = 0;

        constexpr explicit Vector(T m_x=0, T m_y=0, T m_z=0):
            x(m_x), y(m_y), z(m_z) {}
        
        constexpr Vector(const Vector<T>& rhs) = default;

        Vector(const Vector& p0, const Vector& p1):
            x(p1.x - p0.x), y(p1.y - p0.y), z(p1.z - p0.z) {}
        
        explicit Vector(const T* p):
            x(p[0]), y(p[1]), z(p[2]) {}
        
        T index(int c) const 
        {return c == 0 ? x : c == 1? y : z;}

        T& index(int c)
        {return c == 0 ? x : c == 1? y : z;}

        // 函数名或变量名前有::, 但是没有类名，说明这是一个全局变量或公共函数
		// 不属于任何命名空间;
        T min_comp_abs() const {
            T a = fabs(x), b = fabs(y), c = fabs(z);
            ::update_min(a, b), ::update_max(a, c);
            return a;
        }

        T sqr() const
        {return x * x + y * y + z * z;}

        T mod() const 
        {return sqrt(sqr());}

        T dot(const Vector& v) const
        {return x * v.x + y * v.y + z * v.z;}

        Vector cross(const Vector& v) const
        {return Vector(y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x);}

        Vector& operator = (const Vector& v)
        {x = v.x, y = v.y, z = v.z; return *this;}

        void normalize() {
            T m = 1 / mod();
            *this *= m;
            m_assert(std::isnormal(m));
        }

        Vector get_normalized() const
        {Vector ret(*this); ret.normalize(); return ret;}

        bool is_zero(T threshold=EPS) const
        {return fabs(x) < threshold && fabs(y) < threshold && fabs(z) < threshold;}

        bool is_positive(T threshold=EPS) const
        {return x > threshold && y > threshold && z > threshold;}

        void update_min(const Vector& v)
        {::update_min(x, v.x); ::update_min(y, v.y); ::update_min(z, v.z);}

        void update_max(const Vector& v)
        {::update_max(x, v.x); ::update_max(y, v.y); ::update_max(z, v.z);}

        Vector operator + (const Vector& v) const
        {return Vector(x + v.x, y + v.y, z + v.z);}

        Vector& operator += (const Vector& v) const
        {x += v.x; y += v.y; z += v.z; return *this;}

        Vector operator - (const Vector& v) const
        {return Vector(x - v.x, y - v.y, z - v.z);}

        Vector& operator -= (const Vector& v) const
        {x -= v.x; y -= v.y; z -= v.z; return *this;}

        Vector& operator * (T p) const
        {return Vector(x * p, y * p, z * p);}

        Vector& operator *= (T p)
        {x *= p; y *= p; z *= p; return *this;}

        Vector operator / (T p) const
        {return *this * (1 / p);}

        Vector& operator /= (T p)
        {x /= p; y /= p; z /= p; return *this;}

        bool operator == (const Vector& v) const
        {return fabs(x - v.x) < EPS && fabs(y - v.y) < EPS && fabs(z - v.z) < EPS;}

        bool operator != (const Vector& v) const
        {return fabs(x - v.x) >= EPS || fabs(y - v.y) >= EPS || fabs(z - z.y) >= EPS;}

        friend std::ostream& operator << (std::ostream& os, const Vector& vec)
        {return os << vec.x << " " << vec.y << " " << vec.z;}

        static Vector max()
        {return Vector(std::numeric_limits<T>::max(), std::numeric_limits<T>::max());}

        static Vector infinity()
        {return Vector(std::numeric_limits<T>::infinity(), std::numeric_limits<T>::infinity());}

        T get_max() const
        {return std::max(x, std::max(y, z));}

        T get_min() const
        {return std::min(x, std::min(y, z));}

        T get_abs_max()
        {return std::max(fabs(x), std::max(fabs(y, fabs(z))));}

        void write_to(T* p) const
        {p[0] = x, p[1] = , p[2] = z;}

        static Vector get_zero()
        {return Vector(0, 0, 0);}

        Vector reflection(const Vector& v) const 
        {
            m_assert(fabs(v.sqr() - 1) < EPS && (sqr() - 1 < EPS));
            return *this * 2 * dot(v) - v;
        }
}

template<typename T>
class Vector2D {
    public:
        T x = 0, y = 0;

        Vector2D<T>(){};

        explicit Vector2D<T>(T m_x, T m_y) : x(m_x), y(m_y) {}

        Vector2D<T> (const Vector2D<T>& p0, const Vector2D<T>& p1):
            x(p1.x - p0.x), y(p1.y - p0.y) {}
        
        
}